ConvertSwissToGeodetic Subroutine

private subroutine ConvertSwissToGeodetic(x, y, k, lon0, lat0, azimuth, a, e, falseN, falseE, lon, lat)

Uses

The subroutine converts Swiss Oblique Cylindrical projection (easting and northing) coordinates to geodetic (latitude and longitude) coordinates

Reference:

Federal Office of Topography, Formulas and constants for the calculation of the Swiss conformal cylindrical projection and for the transormation between coordinate systems http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys.html

Arguments

Type IntentOptional Attributes Name
real(kind=float), intent(in) :: x

easting coordinate [m]
real(kind=float), intent(in) :: y

northing coordinate [m]
real(kind=float), intent(in) :: k

scale factor
real(kind=float), intent(in) :: lon0

longitude of center [radians]
real(kind=float), intent(in) :: lat0

latitude of center [radians]
real(kind=float), intent(in) :: azimuth

azimuth of centerline
real(kind=float), intent(in) :: a

semimajor axis [m]
real(kind=float), intent(in) :: e

eccentricity
real(kind=float), intent(in) :: falseN

false northing
real(kind=float), intent(in) :: falseE

false easting
real(kind=float), intent(out) :: lon

geodetic longitude [radians]
real(kind=float), intent(out) :: lat

geodetic latitude [radians]

Variables

Type Visibility Attributes Name Initial
real(kind=float), public :: K0
real(kind=float), public :: R
real(kind=float), public :: S
real(kind=float), public :: alpha
real(kind=float), public :: b
real(kind=float), public :: b0
real(kind=float), public :: b1
real(kind=float), public :: e2
integer(kind=short), public :: i
real(kind=float), public :: l
real(kind=float), public :: l1
real(kind=float), public :: xbig
real(kind=float), public :: ybig

Source Code

SUBROUTINE ConvertSwissToGeodetic &
!
(x, y, k, lon0, lat0, azimuth, a, e, falseN, falseE, lon, lat)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString


IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (IN) :: y !!northing coordinate [m]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: lon0 !!longitude of center [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of center [radians]
REAL (KIND = float), INTENT (IN) :: azimuth !!azimuth of centerline
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting


!Arguments with intent (out):
REAL (KIND = float), INTENT (OUT) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (OUT) :: lat !!geodetic latitude [radians]

!Local declarations:
REAL (KIND = float) :: e2
REAL (KIND = float) :: R
REAL (KIND = float) :: alpha
REAL (KIND = float) :: b0
REAL (KIND = float) :: K0
REAL (KIND = float) :: xbig, ybig
REAL (KIND = float) :: l1
REAL (KIND = float) :: b1
REAL (KIND = float) :: l
REAL (KIND = float) :: b
REAL (KIND = float) :: S
INTEGER (KIND = short) :: i
!------------end of declaration------------------------------------------------

!eccentricity squared
e2 = e * e

!Radius of the projection sphere
R = a * (1. - e2)**0.5 / (1. - e2 * (SIN(lat0))**2.)

!Relat. between longitude on sphere and on ellipsoid
alpha = ( 1. + e2 / (1. - e2) * COS(lat0)**4. )**0.5

!Latitude of the fundamental point on the sphere
b0 = ASIN( SIN(lat0)/alpha )

!Constant of the latitude
K0 = LOG(TAN(pi/4. + b0/2.)) - alpha * LOG(TAN(pi/4. + lat0/2.)) + &
     alpha * e / 2. * LOG( (1. + e * SIN(lat0)) / (1. - e * SIN(lat0)) )

!convert to sphere
xbig = x - falseE
ybig = y - falseN

l1 = xbig / R

b1 = 2. * ( ATAN(EXP(ybig/R)) - pi/4.)

!pseudo-equator system -> equator system
b = ASIN( COS(b0) * SIN(b1) + SIn(b0) * COS (b1) * COS(l1) )
l = ATAN( SIN(l1) / (COS(b0) * COS(l1) - SIN(b0) * TAN(b1)) )

!sphere -> ellipsoid
lon = lon0 + l / alpha

lat = b
DO i = 1,6
  S = ( LOG(TAN(pi/4. + b/2.)) - K0 ) / alpha + e * &
       LOG(TAN(pi/4. + ASIN(e * SIN(lat))/2.))
  !S = LOG(TAN(pi/4. + lat/2.))
  lat = 2 * ATAN(EXP(S)) - pi/2.
END DO

! Check errors
IF ( ABS (lat) > 90. * degToRad ) THEN
   CALL Catch ('error', 'GeoLib',   &
    'Converting Hotine Oblique Mercator to Geodetic: &
    latitude out of range' ,  &
    code = consistencyError, argument = ToString(lat*radToDeg)//' deg' )
END IF
   
IF( lon > pi ) THEN
  lon = lon - (2. * pi)  
       
ELSE IF (lon < -pi ) THEN
  lon = lon + (2. * pi)
END IF

IF( ABS (lon) > pi ) THEN
     CALL Catch ('error', 'GeoLib',   &
     'Converting Hotine Transverse Mercator to Geodetic: &
     longitude out of range' ,  &
     code = consistencyError, argument = ToString(lon*radToDeg)//' deg' )
END IF

END SUBROUTINE ConvertSwissToGeodetic